Anima

A C++ library with a set of Python scripts for medical image processing

A. Stamm

Department of Mathematics Jean Leray, UMR CNRS 6629, Nantes University, Ecole Centrale de Nantes, France

2024-11-04

Developing code in Anima

Anima Dependencies1

Overview of tools in Anima

The math-tools module

Your turn

  1. Fork the Anima-Public repository;
  2. Clone your fork of the Anima-Public repository:
# cd somewhere/Anima
git clone https://github.com/astamm/Anima-Public src/
  1. Create a branch for this tutorial:
cd src
git checkout -b anima-tuto
git push --set-upstream origin anima-tuto

My first tool

Getting familiar with data IO

// Content of `math-tools/common_tools/convert_image/animaCreateImage.cxx`

#include <tclap/CmdLine.h>
#include <iostream>
#include <string>

#include <itkImage.h>
#include <itkVectorImage.h>
#include <itkImageRegionIterator.h>

#include <animaReadWriteFunctions.h>

int main(int argc, char **argv)
{
    TCLAP::CmdLine cmd("INRIA / IRISA - VisAGeS/Empenn Team", ' ',ANIMA_VERSION);

    TCLAP::ValueArg<std::string> geomArg("g","geometryfile","Geometry image",true,"","Geometry image",cmd);
    TCLAP::ValueArg<std::string> outArg("o","outputfile","output image",true,"","output image",cmd);
    TCLAP::ValueArg<unsigned int> vdimArg("v","vdim","Force vdim to this value",false,1,"Force vdim",cmd);
    TCLAP::ValueArg<unsigned int> valueArg("b","buffervalue","Value to fill buffer",false,0,"Value to fill buffer",cmd);

    TCLAP::SwitchArg vecArg("V","isvec","Input image is a vector / tensor image",cmd,false);

    try
    {
        cmd.parse(argc,argv);
    }
    catch (TCLAP::ArgException& e)
    {
        std::cerr << "Error: " << e.error() << "for argument " << e.argId() << std::endl;
        return EXIT_FAILURE;
    }

    using ImageType = itk::Image <double,3>;
    using VectorImageType = itk::VectorImage <double,3>;

    bool isVect = vecArg.isSet();
    unsigned int fvdim = vdimArg.getValue();
    if (fvdim > 1)
        isVect = true;

    if (!isVect)
    {
        ImageType::Pointer geomImage = anima::readImage <ImageType> (geomArg.getValue());

        ImageType::Pointer resImage = ImageType::New();
        resImage->Initialize();
        resImage->SetRegions(geomImage->GetLargestPossibleRegion());
        resImage->SetSpacing(geomImage->GetSpacing());
        resImage->SetOrigin(geomImage->GetOrigin());
        resImage->SetDirection(geomImage->GetDirection());

        resImage->Allocate();
        resImage->FillBuffer(valueArg.getValue());

        anima::writeImage <ImageType> (outArg.getValue(),resImage);
    }
    else
    {
        VectorImageType::Pointer geomImage = anima::readImage <VectorImageType> (geomArg.getValue());

        VectorImageType::Pointer resImage = VectorImageType::New();
        resImage->Initialize();
        resImage->SetRegions(geomImage->GetLargestPossibleRegion());
        resImage->SetSpacing(geomImage->GetSpacing());
        resImage->SetOrigin(geomImage->GetOrigin());
        resImage->SetDirection(geomImage->GetDirection());

        unsigned int vdim = geomImage->GetNumberOfComponentsPerPixel();
        if (fvdim > 1)
            vdim = fvdim;

        resImage->SetNumberOfComponentsPerPixel(vdim);
        resImage->Allocate();

        itk::VariableLengthVector <double> tmpData;
        tmpData.SetSize(vdim);
        for (unsigned int i = 0;i < vdim;++i)
            tmpData[i] = valueArg.getValue();

        itk::ImageRegionIterator <VectorImageType> tmpIt(resImage,resImage->GetLargestPossibleRegion());

        while(!tmpIt.IsAtEnd())
        {
            tmpIt.Set(tmpData);
            ++tmpIt;
        }

        anima::writeImage <VectorImageType> (outArg.getValue(),resImage);
    }

    return EXIT_SUCCESS;
}
1
User input arguments are captured by argc and argv and parsed by TCLAP using the CmdLine class. Arguments are declared by type (e.g. ValueArg<double> for a floating-point number, etc.). A variety of other argument types can be used.
2
Use of the TCLAP API requires to include this header.
3
Images are handled by ITK through dedicated classes. These classes are templated on the scalar value type in each voxel and the image dimension. Here we define shortcuts for 3-dimensional scalar and vector images of type double.
4
The itk::Image and itk::VectorImage classes are defined in these headers which therefore have to be included.
5
We read images in and write images out using anima::readImage and anima::writeImage respectively, which are templated on the type of image to handle. The reader expects a string describing the location of the image and outputs a pointer to the memory location where the input image is. The writer expects (i) a string describing the location to which the image should be savec and (ii) the pointer to the memory location where the output image is.
6
We need to include this header to use anima::readImage and anima::writeImage.
7
This code shows how to create an image from scratch using a reference image to copy its geometry. First, the pointer is created and initialized. Then regions, spacing, origin and direction informations are copied from the reference image. Finally, the image is allocated in memory and filled with the desired value.
8
User arguments that have been parsed can be accessed via the .getValue() method for the ValueArg class and via the .isSet() method for the SwitchArg class.
9
To create a vector image, the code is pretty much the same expect that we use the VectorImageType image type.
10
The only subtlety is that now each voxel contains a vector-valued entry which is of type itk::VariableLengthVector. Its length can be retrieved from the reference image via the GetNumberOfComponentsPerPixel().
11
For vector images, there is no FillBuffer() method. Instead, one declares an iterator over the whole region of the image and set its value via the Set() method at each iterator increment until reaching the end of the image region.
12
Iterators over images are ITK classes that are found in specific headers that need to be included.
13
Exit your main program with the EXIT_SUCCESS macro if it runs smoothly up to the end.

Your turn

  1. Copy-paste the previous code in a file that you call animaDoubleLogPValue.cxx;
  2. Modify the code in order to receive as input a 3-dimensional scalar image of p-values and return the image with the double logarithm transform of the p-value in it. We recall that this transformation is

\[ \log \left( - \log(p) \right); \]

  1. Add an optional argument to provide a mask image so that the transformation is applied only on those voxels for which the mask is set to 1.

Tip

Get inspiration from (copy-paste) existing code and modify the minimum to suit your needs!

How to compile this code?

  1. Decide where your tool will end up in the current anima source code tree;
  2. Go there and create a dedicated folder with the name of the tool without the anima prefix and using snake_case formatting (e.g. double_log_pvalue);
  3. At this same location, edit the CMakeLists.txt file by adding add_subdirectory(double_log_pvalue). This tells cmake to add the subdirectory for compilation;
  4. Put the file you wrote (animaDoubleLogPValue.cxx) in the double_log_pvalue folder.

CMakeLists.txt

Now compile your code.

Caution

If you try to compile this code, it will not work because your tool is missing its own CMakeLists.txt file to give instructions to cmake on how to compile our tool (i.e. which libraries should it link with, which headers, etc.).

Do I have to learn cmake as well?

Tip

As usual, let us adopt the minimum effort strategy and see if the existing CMakeLists.txt for animaCreateImage could fit our needs.

Reuse existing code

# Context of `math-tools/common_tools/create_image/CMakeLists.txt`

if(BUILD_TOOLS)

project(animaCreateImage)

## #############################################################################
## List Sources
## #############################################################################

list_source_files(${PROJECT_NAME}
  ${CMAKE_CURRENT_SOURCE_DIR}
  )

## #############################################################################
## add executable
## #############################################################################

add_executable(${PROJECT_NAME}
  ${${PROJECT_NAME}_CFILES}
  )

## #############################################################################
## Link
## #############################################################################

target_link_libraries(${PROJECT_NAME}
  ${ITKIO_LIBRARIES}
  )

## #############################################################################
## install
## #############################################################################

set_exe_install_rules(${PROJECT_NAME})

endif()
1
The if statement checks if the BUILD_TOOLS flag has been set during project configuration via cmake and only compiles this tool if it is the case.
2
This line gives the name of the binary that will be compiled. You need to modify it for your own needs.
3
This part has the necessary code to list sources.
4
This part has the necessary code to add the executable to the list of executables to be compiled.
5
This part has the necessary code to tell the compiler which libraries it should link to. You might need to modify it for your own needs. Check which headers you include in your code to get it right. Most of the time linking with ${ITKIO_LIBRARIES} will be enough.
6
This part has the necessary code to install the binary.
7
Mind the endif() line at the end of the file to close the statement.

Back to your first tool

  1. Copy-paste math-tools/common_tools/create_image/CMakeLists.txt in the double_log_pvalue folder;
  2. Modify its content as needed.
  3. Compile your code!

Note

Anima handles most of the compilation settings for you. You only need to care about the specific CMakeLists.txt for your tool but, even there, most of the time, copy-pasting an already existing one will be enough.

Image Filters

  • A good way of writing clean code by organizing the code that handles the computation in each voxel repeatedly in separate files, out from your main() function;
  • Anima image filters build upon ITK image filters and therefore come with the same benefits such as multi-threading (automatic parallel processing).

From ITK to Anima filters

  • ITK filters transform one or several input images into one or several output images;
  • They feature many advantages such as multi-threading.

Reuse reuse reuse reuse

# Use a terminal to search for existing filters
# based on the one you want to use.

stamm-a@aymerics-mbp Anima % grep -r "animaMasked" *
diffusion/dti/dti_estimator/animaDTIEstimationImageFilter.h:#include <animaMaskedImageToImageFilter.h>
diffusion/dti/flip_tensors/animaFlipTensorImageFilter.h:#include <animaMaskedImageToImageFilter.h>
diffusion/dti/dti_non_central_chi_estimator/animaDTINonCentralChiEstimationImageFilter.h:#include <animaMaskedImageToImageFilter.h>
diffusion/mcm_estimator/animaMCMEstimatorImageFilter.h:#include <animaMaskedImageToImageFilter.h>
math-tools/statistics/animaLocalPatchCovarianceDistanceImageFilter.h:#include <animaMaskedImageToImageFilter.h>
math-tools/statistics/animaLocalPatchMeanDistanceImageFilter.h:#include <animaMaskedImageToImageFilter.h>
math-tools/common/animaMaskedImageToImageFilter.hxx:#include "animaMaskedImageToImageFilter.h"
math-tools/common/animaMaskedImageToImageFilter.h:#include "animaMaskedImageToImageFilter.hxx"
math-tools/statistical_tests/animaCramersTestImageFilter.h:#include <animaMaskedImageToImageFilter.h>
math-tools/statistical_tests/nlmeans_patient_to_group_comparison/animaNLMeansPatientToGroupComparisonImageFilter.h:#include <animaMaskedImageToImageFilter.h>
math-tools/statistical_tests/patient_to_group_comparison/animaPatientToGroupComparisonImageFilter.h:#include <animaMaskedImageToImageFilter.h>
quantitative-mri/estimation/animaT1SERelaxometryEstimationImageFilter.h:#include <animaMaskedImageToImageFilter.h>
quantitative-mri/estimation/gamma_mixture_t2_estimation/animaGammaMixtureT2RelaxometryEstimationImageFilter.h:#include <animaMaskedImageToImageFilter.h>
quantitative-mri/estimation/gmm_t2_estimation/animaGMMT2RelaxometryEstimationImageFilter.h:#include <animaMaskedImageToImageFilter.h>
quantitative-mri/estimation/animaT2EPGRelaxometryEstimationImageFilter.h:#include <animaMaskedImageToImageFilter.h>
quantitative-mri/estimation/combined_relaxometry_estimation/animaCombinedRelaxometryEstimationImageFilter.h:#include <animaMaskedImageToImageFilter.h>
quantitative-mri/estimation/multi_t2_estimation/animaMultiT2RelaxometryEstimationImageFilter.h:#include <animaMaskedImageToImageFilter.h>
quantitative-mri/estimation/animaT1RelaxometryEstimationImageFilter.h:#include <animaMaskedImageToImageFilter.h>
quantitative-mri/estimation/animaT2RelaxometryEstimationImageFilter.h:#include <animaMaskedImageToImageFilter.h>
registration/resamplers/animaOrientedModelBaseResampleImageFilter.h:#include <animaMaskedImageToImageFilter.h>
segmentation/staple/animaMultiThreadedSTAPLEImageFilter.h:#include <animaMaskedImageToImageFilter.h>
segmentation/tissues_em_classification/animaTissuesEMClassificationImageFilter.h:#include <animaMaskedImageToImageFilter.h>

Anatomy of a filter - General

  • An image filter is a C++ class templated over the type of input images and the type of output images;

  • It is made of two files:

    • a .h header file containing method declarations;
    • a .hxx file containing method implementations.

  • Naming convention is animaToolNameImageFilter.h and animaToolNameImageFilter.hxx.

Anatomy of a filter (.h file)

// Content of `math-tools/statistics/animaCramersTestImageFilter.h`

#pragma once

#include <iostream>
#include <animaMaskedImageToImageFilter.h>
#include <itkVectorImage.h>
#include <itkImage.h>
#include <vector>

namespace anima
{

template <class PixelScalarType>
class CramersTestImageFilter :
        public anima::MaskedImageToImageFilter< itk::VectorImage <PixelScalarType, 3> , itk::Image <PixelScalarType, 3> >
{
public:

    /** Standard class typedefs. */
    using Self = CramersTestImageFilter<PixelScalarType>;
    using TInputImage = itk::VectorImage <PixelScalarType, 3>;
    using TOutputImage = itk::Image <PixelScalarType, 3>;
    using Superclass = anima::MaskedImageToImageFilter< TInputImage, TOutputImage >;
    using Pointer = itk::SmartPointer<Self>;
    using ConstPointer = itk::SmartPointer<const Self>;

    /** Method for creation through the object factory. */
    itkNewMacro(Self);

    /** Run-time type information (and related methods) */
    itkTypeMacro(CramersTestImageFilter, MaskedImageToImageFilter);

    /** Image typedef support */
    using InputImageType = TInputImage;
    using OutputImageType = TOutputImage;
    using InputPixelType = typename InputImageType::PixelType;
    using InputImagePointer = typename InputImageType::Pointer;
    using OutputImagePointer = typename OutputImageType::Pointer;

    /** Superclass typedefs. */
    using MaskImageType = typename Superclass::MaskImageType;
    using MaskImagePointer = typename MaskImageType::Pointer;
    using OutputImageRegionType = typename Superclass::OutputImageRegionType;

    /** Set the number of samples to build the distribution */
    itkSetMacro(NbSamples, unsigned long int);
    itkGetMacro(NbSamples, unsigned long int);

    void SetFirstGroupIndexes(std::vector<unsigned int> &group)
    {
        m_FirstGroup.clear();
        m_FirstGroupSize = group.size();

        for (unsigned int i = 0;i < m_FirstGroupSize;++i)
            m_FirstGroup.push_back(group[i]);
    }

    void SetSecondGroupIndexes(std::vector<unsigned int> &group)
    {
        m_SecondGroup.clear();
        m_SecondGroupSize = group.size();

        for (unsigned int i = 0;i < m_SecondGroupSize;++i)
            m_SecondGroup.push_back(group[i]);
    }

    void AddOutlierMask(MaskImageType *mask)
    {
        m_OutlierMasks.push_back(mask);
    }

protected:
    CramersTestImageFilter()
        : Superclass()
    {
        m_NbSamples = 5000;
        m_FirstGroup.clear();
        m_SecondGroup.clear();
        m_FirstGroupSize = 0;
        m_SecondGroupSize = 0;

        m_UseOutlierMasks = false;
        m_OutlierMasks.clear();
    }

    virtual ~CramersTestImageFilter() {}

    void BeforeThreadedGenerateData() ITK_OVERRIDE;
    void DynamicThreadedGenerateData(const OutputImageRegionType &outputRegionForThread) ITK_OVERRIDE;

    //! Bootstrap samples generator for group comparison
    void GenerateBootStrapSamples();

    //! Actually bootstraps a p-value from the computed distance matrix and group samples
    double BootStrap(vnl_matrix <double> &groupDistMatrix, std::vector <double> &inlierWeights);

    //! Computes the Cramers' statistic knowing the distance matrix and the two groups
    double CramerStatistic(vnl_matrix <double> &grpDistMatrix, std::vector <double> &inlierWeights,
                           unsigned int index);

private:
    ITK_DISALLOW_COPY_AND_ASSIGN(CramersTestImageFilter);

    std::vector < unsigned int > m_FirstGroup, m_SecondGroup;
    std::vector < std::vector < unsigned int > > m_SamplesFirstGroup, m_SamplesSecondGroup;

    unsigned long int m_NbSamples;
    unsigned int m_FirstGroupSize, m_SecondGroupSize;

    std::vector <MaskImagePointer> m_OutlierMasks;
    bool m_UseOutlierMasks;
};

} // end namespace anima

#include "animaCramersTestImageFilter.hxx"

Modify public/protected functions & private attributes as needed to clarify your code.

Anatomy of a filter (.hxx file)

// Content of `math-tools/statistics/animaCramersTestImageFilter.hxx`

#pragma once

#include "animaCramersTestImageFilter.h"

#include <itkImageRegionIterator.h>
#include <itkImageRegionConstIterator.h>
#include <itkTimeProbe.h>
#include <itkProgressReporter.h>

#include <cmath>
#include <ctime>
#include <random>
namespace anima
{

template <class PixelScalarType>
void
CramersTestImageFilter<PixelScalarType>
::BeforeThreadedGenerateData ()
{
    Superclass::BeforeThreadedGenerateData();

    this->GetOutput()->FillBuffer(0);

    // Checking consistency of the data and parameters

    unsigned int nbInputs = this->GetNumberOfIndexedInputs();
    if (nbInputs <= 1)
        itkExceptionMacro("Error: Not enough inputs available... Exiting...");

    if ((m_FirstGroupSize + m_SecondGroupSize) != nbInputs)
        itkExceptionMacro("Groups data not clearly wrong... Exiting...");

    m_SamplesFirstGroup.clear();
    m_SamplesSecondGroup.clear();

    this->GenerateBootStrapSamples();

    m_UseOutlierMasks = (m_OutlierMasks.size() == nbInputs);
}

template <class PixelScalarType>
void
CramersTestImageFilter<PixelScalarType>
::DynamicThreadedGenerateData(const OutputImageRegionType &outputRegionForThread)
{
    using InIteratorType = itk::ImageRegionConstIterator < TInputImage >;
    using OutRegionIteratorType = itk::ImageRegionIterator < OutputImageType >;

    using MaskRegionIteratorType = itk::ImageRegionIterator < MaskImageType >;

    OutRegionIteratorType outIterator(this->GetOutput(), outputRegionForThread);
    MaskRegionIteratorType maskIterator (this->GetComputationMask(), outputRegionForThread);

    std::vector <InIteratorType> inIterators(this->GetNumberOfIndexedInputs());
    for (unsigned int i = 0;i < this->GetNumberOfIndexedInputs();++i)
        inIterators[i] = InIteratorType(this->GetInput(i), outputRegionForThread);

    std::vector <MaskRegionIteratorType> outlierMasksIterators(m_OutlierMasks.size());

    if (m_UseOutlierMasks)
    {
        for (unsigned int i = 0;i < this->GetNumberOfIndexedInputs();++i)
            outlierMasksIterators[i] = MaskRegionIteratorType(m_OutlierMasks[i],outputRegionForThread);
    }

    std::vector < InputPixelType > firstGroupData(m_FirstGroupSize), secondGroupData(m_SecondGroupSize);
    vnl_matrix <double> cramerDistMatrix(this->GetNumberOfIndexedInputs(),this->GetNumberOfIndexedInputs(),0.0);
    std::vector <double> inlierProbabilities(this->GetNumberOfIndexedInputs(),1);
    unsigned int vectorSize = this->GetInput(0)->GetNumberOfComponentsPerPixel();

    while (!outIterator.IsAtEnd())
    {
        if (maskIterator.Get() == 0)
        {
            outIterator.Set(0.0);
            ++outIterator;
            ++maskIterator;

            for (unsigned int i = 0;i < this->GetNumberOfIndexedInputs();++i)
                ++inIterators[i];

            if (m_UseOutlierMasks)
            {
                for (unsigned int i = 0;i < this->GetNumberOfIndexedInputs();++i)
                    ++outlierMasksIterators[i];
            }

            continue;
        }

        // Voxel is in mask, so now let's go for it and compute the stats
        for (unsigned int i = 0; i < m_FirstGroupSize;++i)
        {
            unsigned int indFirst = m_FirstGroup[i];
            firstGroupData[i] = inIterators[indFirst].Get();
        }

        for (unsigned int i = 0; i < m_SecondGroupSize;++i)
        {
            unsigned int indSec = m_SecondGroup[i];
            secondGroupData[i] = inIterators[indSec].Get();
        }

        // Now that data is grouped, compute the distance matrix
        for (unsigned int i = 0; i < m_FirstGroupSize;++i)
        {
            for (unsigned int l = 0;l < m_SecondGroupSize;++l)
            {
                double dist = 0;
                for (unsigned int j = 0;j < vectorSize;++j)
                    dist += (firstGroupData[i][j] - secondGroupData[l][j]) * (firstGroupData[i][j] - secondGroupData[l][j]);

                cramerDistMatrix(i,m_FirstGroupSize + l) = sqrt(dist);
                cramerDistMatrix(m_FirstGroupSize + l,i) = cramerDistMatrix(i,m_FirstGroupSize + l);
            }

            for (unsigned int l = i+1;l < m_FirstGroupSize;++l)
            {
                double dist = 0;
                for (unsigned int j = 0;j < vectorSize;++j)
                    dist += (firstGroupData[i][j] - firstGroupData[l][j]) * (firstGroupData[i][j] - firstGroupData[l][j]);

                cramerDistMatrix(i,l) = sqrt(dist);
                cramerDistMatrix(l,i) = cramerDistMatrix(i,l);
            }
        }

        for (unsigned int i = 0; i < m_SecondGroupSize;++i)
        {
            for (unsigned int l = i+1;l < m_SecondGroupSize;++l)
            {
                double dist = 0;
                for (unsigned int j = 0;j < vectorSize;++j)
                    dist += (secondGroupData[i][j] - secondGroupData[l][j]) * (secondGroupData[i][j] - secondGroupData[l][j]);

                cramerDistMatrix(m_FirstGroupSize + i,m_FirstGroupSize + l) = sqrt(dist);
                cramerDistMatrix(m_FirstGroupSize + l,m_FirstGroupSize + i) = cramerDistMatrix(m_FirstGroupSize + i,m_FirstGroupSize + l);
            }
        }

        if (m_UseOutlierMasks)
        {
            for (unsigned int i = 0;i < this->GetNumberOfIndexedInputs();++i)
                inlierProbabilities[i] = (outlierMasksIterators[i].Get() == 0);
        }

        outIterator.Set(this->BootStrap(cramerDistMatrix,inlierProbabilities));

        ++outIterator;
        ++maskIterator;

        for (unsigned int i = 0;i < this->GetNumberOfIndexedInputs();++i)
            ++inIterators[i];

        if (m_UseOutlierMasks)
        {
            for (unsigned int i = 0;i < this->GetNumberOfIndexedInputs();++i)
                ++outlierMasksIterators[i];
        }

        this->IncrementNumberOfProcessedPoints();
    }
}

template <class PixelScalarType>
void
CramersTestImageFilter<PixelScalarType>
::GenerateBootStrapSamples()
{
    m_SamplesFirstGroup.resize(m_NbSamples + 1);
    m_SamplesSecondGroup.resize(m_NbSamples + 1);

    m_SamplesFirstGroup[0] = m_FirstGroup;
    m_SamplesSecondGroup[0] = m_SecondGroup;

    unsigned int nbFirstGroup = m_FirstGroup.size();
    unsigned int nbSecondGroup = m_SecondGroup.size();
    unsigned int nbSubjects = nbFirstGroup + nbSecondGroup;

    unsigned int minNbGroup = std::min(nbFirstGroup,nbSecondGroup);
    bool isFirstGroupMin = (nbFirstGroup <= nbSecondGroup);

    std::vector <unsigned int> sampleGen, sampleGenOtherGroup;
    std::mt19937 generator(time(0));
    std::uniform_real_distribution<double> unifDistr(0.0, 1.0);

    for (unsigned int i = 1;i <= m_NbSamples;++i)
    {
        sampleGen.clear();
        sampleGenOtherGroup.clear();

        unsigned int j = 0;
        while (j < minNbGroup)
        {
            int tmpVal = std::min((int)(nbSubjects - 1),(int)floor(unifDistr(generator) * nbSubjects));
            if (tmpVal < 0)
                tmpVal = 0;

            bool isAlreadyIndexed = false;
            for (unsigned int k = 0; k < j;++k)
            {
                if (tmpVal == sampleGen[k])
                {
                    isAlreadyIndexed = true;
                    break;
                }
            }

            if (!isAlreadyIndexed)
            {
                sampleGen.push_back(tmpVal);
                j++;
            }
        }

        // New sample gen that is not already here... Add it to the result vectors
        for (unsigned int j = 0;j < nbSubjects;++j)
        {
            bool isAlreadyIndexed = false;
            for (unsigned int k = 0;k < minNbGroup;++k)
            {
                if (j == sampleGen[k])
                {
                    isAlreadyIndexed = true;
                    break;
                }
            }

            if (!isAlreadyIndexed)
                sampleGenOtherGroup.push_back(j);
        }

        if (isFirstGroupMin)
        {
            m_SamplesFirstGroup[i] = sampleGen;
            m_SamplesSecondGroup[i] = sampleGenOtherGroup;
        }
        else
        {
            m_SamplesFirstGroup[i] = sampleGenOtherGroup;
            m_SamplesSecondGroup[i] = sampleGen;
        }
    }
}

template <class PixelScalarType>
double
CramersTestImageFilter<PixelScalarType>
::BootStrap(vnl_matrix <double> &groupDistMatrix, std::vector <double> &inlierWeights)
{
    // First index is the true group separation (see GenerateBootStrap)
    double dataVal = this->CramerStatistic(groupDistMatrix,inlierWeights,0);

    std::vector <double> statsValues(m_NbSamples,0);

    for (unsigned long int i = 0;i < m_NbSamples;++i)
        statsValues[i] = this->CramerStatistic(groupDistMatrix,inlierWeights,i+1);

    std::sort(statsValues.begin(),statsValues.end());

    if (dataVal >= statsValues[statsValues.size() - 1])
        return 0;

    unsigned int position = 0;

    while(position < statsValues.size())
    {
        if (dataVal <= statsValues[position])
            break;

        ++position;
    }

    double resVal = 0;

    if (position > 0)
    {
        resVal = m_NbSamples - position + (statsValues[position - 1] - dataVal) / (statsValues[position] - statsValues[position - 1]);
        resVal /= m_NbSamples;
    }
    else
    {
        // Here statsValues[position - 1] doesn't exist, replacing with 0
        resVal = m_NbSamples - position - dataVal / statsValues[position];
        resVal /= m_NbSamples;
    }

    return resVal;
}

template <class PixelScalarType>
double
CramersTestImageFilter<PixelScalarType>
::CramerStatistic(vnl_matrix <double> &grpDistMatrix, std::vector <double> &inlierWeights,
                  unsigned int index)
{
    unsigned int nbFirstGroup = m_SamplesFirstGroup[index].size();
    unsigned int nbSecondGroup = m_SamplesSecondGroup[index].size();

    double firstTerm = 0, secondTerm = 0, thirdTerm = 0;
    double Z1 = 0, Z2 = 0, Z3 = 0;

    for (unsigned int i = 0; i < nbFirstGroup;++i)
    {
        unsigned int indFirst = m_SamplesFirstGroup[index][i];
        for (unsigned int j = 0;j < nbSecondGroup;++j)
        {
            double alpha = inlierWeights[indFirst]*inlierWeights[m_SamplesSecondGroup[index][j]];
            firstTerm += alpha*grpDistMatrix(indFirst,m_SamplesSecondGroup[index][j]);
            Z1 += alpha;
        }

        for (unsigned int j = i+1;j < nbFirstGroup;++j)
        {
            double alpha = 0.5*inlierWeights[indFirst]*inlierWeights[m_SamplesFirstGroup[index][j]];
            secondTerm += 2*alpha*grpDistMatrix(indFirst,m_SamplesFirstGroup[index][j]);
            Z2 += 2*alpha;
        }

        // Don't forget to add the diagonal term for Z2
        Z2 += inlierWeights[indFirst]*inlierWeights[m_SamplesFirstGroup[index][i]];
    }

    firstTerm /= Z1;
    secondTerm /= (2*Z2);

    for (unsigned int i = 0; i < nbSecondGroup;++i)
    {
        int indSec = m_SamplesSecondGroup[index][i];
        for (unsigned int j = i+1;j < nbSecondGroup;++j)
        {
            double alpha = inlierWeights[indSec]*inlierWeights[m_SamplesSecondGroup[index][j]];
            thirdTerm += 2*alpha*grpDistMatrix(indSec,m_SamplesSecondGroup[index][j]);
            Z3 += 2*alpha;
        }

        // Don't forget to add the diagonal term for Z3
        Z3 += inlierWeights[indSec]*inlierWeights[m_SamplesSecondGroup[index][i]];
    }

    thirdTerm /= (2*Z3);

    return firstTerm - secondTerm - thirdTerm;
}

} // end namespace anima
1
This method is used to perform any computation that can be done once to pre-compute values that can be later used for computations in each voxel.
2
This method is used by each worker to process the voxels that have been assigned to it. It receives the appropriate region for this purpose. Iterators are therfore instantiated within this method to cover the appropriate region.

Back again to our first tool

  1. Write an image filter to handle the p-value transformation;

    • Write animaDoubleLogPValueImageFilter.h;
    • Write animaDoubleLogPValueImageFilter.hxx;

  2. Modify animaDoubleLogPValue.cxx to call the filter instead of doing the computation directly in the main() function;

  3. Add an optional argument to the parser for setting the number of cores to parallelize over.

Tip

Again, get inspiration from the existing code to understand how to modify the main() function. See for instance animaCramersTest.cxx.